Event Horizons in Numerical Relativity II: Analyzing the Horizon 
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We present techniques and methods for analyzing the dynamics of event horizons in numerically 
constructed spacetimes. There are three classes of analytical tools we have investigated. The 
first class consists of proper geometrical measures of the horizon which allow us comparison with 
perturbation theory and powerful global theorems. The second class involves the location and 
study of horizon generators. The third class includes the induced horizon 2-metric in the generator 
comoving coordinates and a set of membrane-paradigm like quantities. Applications to several 
distorted, rotating, and colliding black hole spacetimes are provided as examples of these techniques. 



I. INTRODUCTION 

Black holes play an important role in general relativity 
and astrophysics. They are characterized both by space- 
time singularities within them and by their horizons that 
cover the singularities from the outside world. In this pa- 
per we develop a set of tools for analyzing the dynamics 
of black hole horizons. 

The event horizon (EH) of a black hole is defined as the 
boundary of the causal past of future null infinity X + . As 
such the EH surface is traced out by light rays that never 
reach future null infinity and never fall into the black 
hole singularity. This surface responds to infalling matter 
and radiation and to the gravitational fields of external 
bodies. In the membrane paradigm of black holes, the 
horizon fully characterizes the dynamical interactions of 
a black hole with its surroundings Q. The important 
role of the horizon in the study of black holes motivates 
us to carry out a systematic study of horizon dynamics 
in numerical relativity. 

While much work has been done on the properties 
of stationary black holes and small perturbations about 
them, little is known about the properties of highly dy- 
namical black hole spacetimes. For example, the cos- 
mic censorship conjecture j|, which suggests that space- 
time singularities should be clothed by event horizons, 
demands study into the existence of horizons. The hoop 
conjecture ^Qj, which states that a black hole horizon 
forms if and only if a matter source becomes sufficiently 
compact in all directions, begs the question of how spher- 
ical must a black hole horizon be. Caustics, or singular 
points in the congruence of photons tracing out the hori- 
zon where new generators can join the horizon, can occur, 
but under what conditions do they appear? And what 
are the properties of these caustics? One would also like 
to know to what extent one can understand interactions 



of black holes with their astrophysical environment in 
terms of properties of the EH. Studies of most of these 
questions have to date only been made in very idealized 
circumstances or quasi-stationary spacetimes. But as- 
pects of each of these open questions are amenable to 
study with the numerical methods we describe. 

Due to strong field nonlinearities, black hole horizons 
are difficult to study analytically. Therefore we turn to 
numerical treatment which is now routinely able to gen- 
erate highly dynamical, axisymmetric black hole space- 
times evolved beyond t = lOOAf , where M is the ADM 
mass of the spacetime. Many such axisymmetric studies 
of highly distorted rotating and non-rotating black holes 
and colliding black holes have been performed in recent 
years || |^| . Three dimensional black hole evolutions are 
approaching the accuracy of axisymmetric calculations 
[p|-|T3||. Together with the ability to find and analyze 
event horizons, these simulations provide us with a new 
opportunity to study black hole dynamics. 

We recently proposed methods for the study of the EH 
in numerically generated spacetimes In a series of 

followup papers, we give details of the methods and their 
applications to various black hole spacetimes. The first 
paper in this series referred to hereafter as Paper I, 
detailed the method for locating the EH in a dynamical 
spacetime, and showed the high degree of accuracy with 
which the EH can be located. In this second paper, we 
focus on the tools constructed for analyzing the dynamics 
of the EH. 

There are several aims of the present paper. We show 
three different sets of tools that can be used to analyz- 
ing the dynamics of the EH, and how one can construct 
them in numerical relativity. We show how accurately 
the quantities used in these tools can be constructed with 
present numerically generated black hole spacetimes. We 
demonstrate the applicability of these tools to various 
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spacetimes of interest. In fact, these tools apply imme- 
diately to almost all numerically generated black hole 
spacetimes we have constructed to date. This paper de- 
scribes the tools that elucidate the physics of the EH, 
and the accuracy with which we can (or cannot) evalu- 
ate these measures; the emphasis is not on the physics 
itself. The physics we learn using these tools will be dis- 
cussed in a later paper in this series. 

We have developed and present three sets of tools for 
analyzing the EH. First, we present a set of geometric 
measures of the horizon as a two dimensional surface 
in a curved 3D space-like slice of constant time. These 
tools include proper circumferences, proper area, Gaus- 
sian curvature, the embedding of the surface in Euclidean 
space, and the embedding history. Second, we discuss 
how the horizon generators can be constructed. This 
construction also gives the locus of generators that will 
join the horizon in the future at caustic points on the 
horizon surface. Third, we present a set of tools from 
the membrane paradigm of black holes Q for analyzing 
the generators and the physics they contain, such as the 
horizon 2-metric in generator co-moving coordinates, 7^ 
and quantities derived from and connected to it, such as 
the expansion <d H , shear a^ b , surface gravity gn : and 
Hajicek field fif . 

To illustrate the use of these horizon tools, we ap- 
ply them to several spacetimes. We consider the 
Schwarzschild and Kerr analytic black hole spacetimes 
to show the basic principles involved, and to test the 
accuracy of the methods. Also, we apply them to fully 
nonlinear, highly dynamical black hole systems, such as a 
distorted Schwarzschild BH, a distorted Kerr BH, and the 
collision of two black holes (the Misner data [^6|). Our 
tools can be applied to almost all numerical black hole 
spacetimes we have presently constructed, and should be 
applicable to future black hole spacetimes as well. 

The structure of this paper is as follows: In Sec. || we 
briefly review the method we developed to find the loca- 
tion of the EH. In Sec. [II we show various ways to ex- 
tract important information from the EH surface location 
in the spacetime, including studying the topology, area, 
various circumferences, Gaussian curvature, and geomet- 
ric embeddings of the surface. In Sec. IV we show how to 
find the actual generators of the EH, and the information 
their paths can bring. In Sec. [v| we discuss how one can 
apply ideas developed in the membrane paradigm Q to 
numerically generated black hole spacetimes. Through- 
out the paper, we illustrate these ideas with examples 
from numerically generated black hole spacetimes. 



II. LOCATING THE EH IN A NUMERICALLY 
GENERATED BLACK HOLE SPACETIME 

Our method for locating event horizons in numerical 
relativity was detailed in Paper I. In order to define our 
notation, and because our analysis here is closely related 



to our EH finding method, we briefly review it here. The 
essence of the EH finding method can be summarized in 
four steps: 

(i) At late times after the dynamical evolution we seek 
to analyze (that is, when the black hole spacetime has 
returned to approximate stationarity; e.g., after the coa- 
lescence of two black holes or after all incident gravitation 
radiation has either radiated into the hole or into the far 
wave zone) the position of the EH can often be located 
approximately. We can identify a region of the late-time 
spacetime which contains the EH, which we call the hori- 
zon containing domain (HCD). 

(ii) We trace the evolution of the HCD backward in time 
by tracing its outer and inner boundaries as null surfaces. 
A function describing a null surface at t = tj, x l = x l j 



f(t = t f ,x})=0, 
satisfies the equation 

d,fd»f = 0, 



(1) 
(2) 



or, 



a , -g u dif + VOW) 2 - y'v^„/vv 

d *f = Jt • ( 3 ) 

for outgoing surfaces. The fact that this method repre- 
sents the position of the EH directly as a function f(t, x l ) 
is particularly convenient in our construction of horizon 
analysis tools, as we shall see below. 
(m)The strength of our method stems from the fact that 
the inner and outer boundaries of the HCD converge 
together quickly when integrated backwards in time in 
many cases of interest. When the distance between the 
two boundaries in a time slice becomes significantly less 
than the grid separation used in the construction of the 
spacetime, we have accurately located the EH. This con- 
dition can often be met through the entire regime of in- 
terest here. 

(iv) The choice of parameterization of the surface is im- 
portant. For the axisymmetric spacetimes used as exam- 
ples in this paper, one convenient choice is 



/ = T) - S(t, ( 



(4) 



where 77 is a radial coordinate and 9 is polar angular 
coordinate. In what follows, we assume that this func- 
tion f(t, 77, 6) has been obtained for the numerically con- 
structed spacetime. 

In axisymmetric two black hole spacetimes, such as 
those generated in j7 17, we can use the parameteri- 
zation in Eq. (0) to trace the event horizon through the 
merger phase. In the Cadez coordinate system |lljj where 
coordinates are centered around each individual throat 
and the axis below the throat is a line of constant 9, this 
parameterization allows us to trace a single hole by ap- 
plying an upwinded condition on the horizon at the axis 
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before coalescence. In the recently proposed "Class I" co- 
ordinate system fL8f where the coordinates are centered 
around the throat and the axis, forming peanut shaped 
radial coordinate lines near the throats, this parameter- 
ization will represent a null surface which contains both 
the horizon and the null surface which represents the lo- 
cus of generators waiting to join the horizon; a simple 
symmetry boundary condition on the equator suffices. 
The locus can also be located in the Cadez system by 
using an alternate (p, z) parameterization, as described 
in Paper I. We will use simulations generated in both 
coordinate systems interchangeably here. 



III. STUDYING THE EH SURFACE 

A. Topology of the EH 

We note that one important, but easy to obtain piece 
of information contained in f(t, x 1 ) — is the topology 
of the EH at a constant time slice. In this section we 
show an interesting case of the EH undergoing a change 
of topology. We apply our EH finding method above to a 
numerically generated spacetime representing the head- 
on collision of two equal mass black holes with axisym- 
metry. In Fig. [I] we show the function f(t,x l ) = at 
two times for the case of Misner time symmetric initial 
data, described by two throats connecting two identical 
asymptotically flat sheets, evolved by a code described 
in Ref. |?J (the "Cadez" code). The case considered here 
is for the Misner parameter /i = 2.2, for which the initial 
distance between the throats is 8.92M, where M is the 
1/2 the ADM mass. For details of the initial data set, 
see Ref. MM. 
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FIG. 1. We show the topology of the EH for the collision 
of two black holes using data generated with the Cadez code 
from the Misner y, = 2.2 initial data. The solid lines show 
the EH at t = OM and t = 7.5M. The dotted line at t = OM 
shows the locus of generators which will join the horizon in 
the future. We note that between t = OM and t = 7.5M the 
horizon undergoes a non-trivial topology change. We note 
that "crossover" type caustics form on the axis. We shall see 
below that new generators join the horizon at these points. 

At t = the EH has the topology of two disconnected 
two spheres represented by the solid lines centered near 
z = ±1 in the p — z plane. We note that the function 
f(t = 0, x l ) gives not just the location of the EH, but also 
the locus of the future horizon generators before they join 
the horizon. At t = 7.5M the horizon has the topology 
of a single sphere. 

We treat this change of EH topology by following the 
surface function backwards in time. We can trace the 
horizon from i ~ 75 or 100AT to a "dumbell" shaped hori- 
zon at t = 7.5M. In Fig. [[], we start with this "dumbell" 
shaped horizon. Tracing this surface backward towards 
t = OM, we see that the central part of the surface shrinks 
rapidly, and the left and right hand sides cross, indicat- 
ing the change in topology. At t — the portion of the 
surface corresponding to the locus of photons which will 
join the horizon, but have not yet done so, is given as 
a dashed line. As discussed in Sec. IV, the crossing of 
the surface signals that photons are leaving the horizon, 
going backwards in time. The crossed portion of the sur- 
face (shown as a dashed line in Fig|y) is no longer on 
the EH, but represents the surface of horizon generators 
"waiting to be born" , as they will join the horizon at a 
future time. For the work here, we define a point on the 
horizon where generators cross as a caustic, and there- 
fore, at the point where the generators cross and join the 
horizon, the horizon has a caustic point. This caustic at 
the cusp in the event horizon is discussed further below, 
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and also in Ref. po| . 



deviates significantly from the analytic value of 16irM 2 
as the evolution continues. 



B. Geometry of the EH Surface 

The function for the surface f(t, x l ), together with the 
metric induced on the surface, gives the intrinsic geome- 
try of the EH, from which important physical properties 
can be determined. In this section we present a set of 
tools which allow one to study the intrinsic properties of 
the surface. 



1. Area 

There has been extensive study of the surface area of 
black hole event horizons in general relativity pi| . The 
area plays a central role in the thermodynamics of black 
holes. As area is a quantity directly used in analytic stud- 
ies, it is important to be able to study the dynamical evo- 
lution of the area of the EH in a numerically constructed 
spacetime, both for understanding the spacetime, and 
also as a diagnostic tool for the accuracy of the numeri- 
cal treatment. 

Construction of the surface area as function of time 
is straightforward. Here we show how one computes the 
area mainly for establishing the notation used in this pa- 
per. A surface f(to,x l ) = determines the coordinate 
location x l = x l (6, 4>), (i — 1, 2, 3), of the surface at time 
to, where we regard the surface as being parametrized by 
two surface coordinates x a — (8,(f>),(a — 1,2). Denote 
the spatial line element of the spacelike hypersurface at 
t = t by 



da' 



9ij 



dx l dx 1 , 



(5) 



and so we can define an induced horizon 2-metric as 



Tab 9ij 



dx l dxi 
dx a dx b 



The surface area at a time t is then given by 



A(t) 



J^dx*, 



(G) 



(7) 



where 7 is the determinant of 7 a {,. 

In Fig. |^ we show A(t) for a Schwarzschild black hole 
evolved with maximal slicing. The dotted line (labeled 
Standard ADM) shows the results obtained using the 
standard numerical treatment as described in Ref. p2| . 
In this case, the calculation is carried out using a ID code 
with 200 grid zones. It is well known that when evolved 
with such slicings and without shift, the EH will expand 
outward in the radial direction in coordinate space. At 
the same time, a sharp peak in the radial metric function 
develops near the EH. Due to numerical error caused by 
the inability to resolve this sharp peak, the function A(t) 
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FIG. 2. We show the area of the EH traced out for a 
Schwarzschild spacetime evolved with two different methods. 
The dotted line shows the EH evolution with a standard 3+1 
ADM evolution scheme. The dashed line shows the evolution 
with an apparent horizon boundary condition. The solid line 
shows the analytic value. We note that the same method is 
used for both spacetimes. The error in the standard ADM 
spacetime is due to inaccuracies in the spacetime metric, not 
the EH finder. 

We compare this result to the case of the dashed line 
in Fig. ^, which is obtained by applying Eq. (0) to a 
Schwarzschild spacetime constructed with the same grid 
parameters, but with an apparent horizon boundary con- 
dition |23|,Q. The improvement in accuracy is dramatic. 

We stress that the issue here is the accuracy of the nu- 
merically constructed spacetime, and not the accuracy 
of the EH finding method; an identical finder is used for 
both the dashed and dotted line. The error in A(t) in the 
case of the dotted line is dominated by the error in the 
spacetime data. The relative numerical error in finding 
the function f(t,r,9,(f>) — as the position of the EH 
is small compared to the errors in the background space- 
time. The agreement of the dashed line with the analytic 
value A(t) — IQttM 2 suggests that the error in A(t) de- 
termined with the apparent horizon boundary condition 
spacetime is only about 1% at t = 100M. Though sim- 
ple, this is an illustrative example of using the horizon 
analysis as tool to understand the accuracy of a given 
numerical spacetime. 



2. Circumference 

For black holes with symmetries, the definitions of 
some circumferences are geometrically meaningful. For 
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example, in axisymmetric spacetimes one can define a 
polar circumference C p , and for spacetimes with a reflec- 
tion symmetry around the equatorial plane, an equato- 
rial circumference C e . In the axisymmetric system, with 
d/d(f> being the azimuthal Killing vector, we can take the 
horizon coordinates x a to be those tied to the symmetry 
axis, 

x a = (M) = (0,0). (8) 

The polar circumference C p , the circumference of a line 
with 4> — const, is 




The equatorial circumference C e , which is the loop 
around the horizon at 9 = tt/2, is given by 

C e = [ ^/ lab dx a dx b (10) 

J6=n/2 

What is often more interesting is not C p or C e by them- 
selves, but their ratio C r = jf-. This defines an effec- 
tive shape parameter for axisymmetric surfaces. Roughly 
speaking, if C r > 1 or C r < 1 the surface is prolate or 
oblate, respectively. 

a. Shape of the Analytic Kerr Horizon In Fig. || we 
show the quantity C r for Kerr black holes with vari- 
ous rotation parameters a. The numerical simulation of 
such spacetimes have been discussed in Refs. p5|,p|,p6[ . 
The Kerr spacetimes we consider here, however, are not 
evolved, but rather the analytic (stationary) Kerr solu- 
tion in the logarithmic radial rj,9 coordinates. The use 
of analytic data enables us to test directly the accuracy 
of our horizon treatment, without being affected by the 
error of representing the spacetime on a numerical grid 
with finite resolution. The solid line shows the analytic 
value , and the diamonds are data points obtained by 
applying our methods to Kerr spacetimes and measuring 
their circumferences as described above. The agreement 
on the plot is excellent. 
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FIG. 3. We plot the ratio CV of the polar to equatorial 
circumferences of the event horizons for Kerr black holes 
with various rotation parameters a. The diamonds show 
data points obtained by applying our method to various Kerr 
spacetimes, while the solid line shows the analytic result. The 
agreement is excellent. 

b. Shape of distorted EHs Next we consider the event 
horizons of highly distorted black holes. An important 
open question about the nature of black holes is the fol- 
lowing: "Is the event horizon always rather spherical, 
as suggested by the Hoop Conjecture |,||,|9|?" This 
question can be addressed to some extent by studying 
the EH of black holes distorted by axisymmetric gravita- 
tional waves. The initial data construction has been de- 
scribed in detail in Ref. |3(J. For our purposes it suffices 
to note that the system corresponds to a time symmet- 
ric torus of gravitational waves, whose amplitude and 
shape are specifiable as parameters, that surround an 
Einstein-Rosen bridge. In Fig. |] we survey the event 
horizons at the initial time t = for a range of black 
hole data sets with fixed Brill wave shape parameters 
((a = 1.0, r/ = 0.0, n = 2) in the language of Ref. [|o)) 
representing a quadrupolar wave centered on the black 
hole throat with a width of order 1M. To find the EH 
at the initial time, we first evolve the initial data to a 
late time, and then trace the EH backwards through the 
evolved data, as described above. Fig. [| shows the EH 
parameter C r for the initial data as a function of the 
Brill wave amplitude Qo. We see that in the range of 
parameters investigated C r can be rather large (almost 3 
in Fig. but does seem to have a maximum in Qq space 
when measured at t = 0. In contrast, at the same inci- 
dent wave amplitude, the AH has much larger amplitude 
and is increasing in increasing amplitude to substantially 
larger distortions than the EH. As this paper is restricted 
to the introduction of the analysis tools of the EH, we will 
defer an exhaustive parameter search and discussion of 
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the physical implications of this result to a later paper, 
including a comparison with Fig. 4 of |^| , where the up- 
per bound on the distortion of the apparent horizon is 
found to be orders of magnitude larger than that of the 
EH. In the initial data, the AH is far inside the EH, but 
after a short evolution in these spacetimes, the AH will 
quickly pop out (the AH is generically spacelike) to be 
closer to (but still inside) the EH. 
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FIG. 4. We plot the shape parameter CV for a series of 
event horizons of non-rotating black holes surrounded by grav- 
itational waves of varying strengths, denoted as triangles. For 
comparison, we show the behavior of the apparent horizon, de- 
noted as diamonds. We note that CV of the AH at Qo — 1.8 
is around 17, in contrast with the EH, which has CV around 
2.5. 

While Fig. || shows C r for highly distorted black hole 
event horizons at t = 0, the same function also provides 
important insight into the evolution of these horizons. In 
particular, we study the case of Q = 1.0. As shown in 
Fig. [|, when this black hole evolves, its horizon evolves 
towards sphericity, overshoots, and oscillates about its 
equilibrium, spherical configuration. The frequency and 
decay rate are, to very high accuracy, the quasinormal 
mode (QNM) of the black hole as determined by pertur- 
bation theory. For comparison, a fit to the two lowest 
QNM frequencies is given by the dashed line. 



o 
ii 

h 

o 



1.04 r 



1.02 



1.00 



0.98 U 

10 20 30 40 50 60 70 




Event Horizon 
QNM Fit 



10 



20 



30 40 

t/M , 



50 



60 



70 



FIG. 5. The time development of the shape parameter CV 
is shown for a single black hole distorted by a Brill wave. We 
note that although the black hole is initially very distorted 
(Cr = 2.9) it quickly settles down to the quasi-normal mode 
ringing predicted by perturbation theory. 

The oscillations of the EH are a common dynamical 
feature of black holes. In Fig. |^ we show the oscillation 
of the EH in the two black hole collision simulated using 
the "Class I" coordinates with the Misner separation pa- 
rameter fi — 2.2, whose coordinate location was shown 
in Fig. |l|. We show the oscillation of the single horizon 
which forms after the merger of the two individual black 
holes. The results are similar to those of the single dis- 
torted black hole described above, showing the generic 
nature of this phenomenon. 
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FIG. 6. The time development of the shape parameter CV 
is shown for two colliding black holes. We note again that 
despite the violent initial beginnings, the system settles down 
to ringing behavior at late times. Note that the definition of 
the shape parameter is not appropriate until after the coales- 
cence. 



3. Gaussian Curvature 

The Gaussian curvature is a local property describing 
the two principal radii of curvature at each point on a 
surface. It has been found to be very useful in analyzing 
the dynamics of the AH and it applies equally well 
to the EH surface. The general formula for the Gaussian 
curvature of a 2-surface in a spacelike hypersurface with 
3-metric is n — 2R where R is the Ricci scalar of the 
2-sphere with the induced 2-metric. 

Fig. |?] shows the time evolution history of the Gaus- 
sian curvature for the highly distorted hole studied in Fig. 
^. Horizon history diagrams like this have proved very 
useful in showing the development of apparent horizon 
surfaces in time (3l]] , and here we apply them to the EH 
surface. The figure shows the evolution of the Gaussian 
curvature as a gray-scale across the surface as a function 
of time (horizontal axis). We use the z— axis embedding 
of the horizon (described below) as the vertical axis in the 
plots, k is larger initially near the equator and then oscil- 
lates between the poles and equator. The checkerboard 
pattern is typical of a predominantly 1 — 2 distortion of 
the horizon, as discussed in Ref. |31|] (there in the AH 
case) . The frequency of oscillation of the horizon surface 
can be read off directly from the figure. We see that it 
has a period of about 17M, which is the fundamental 
period of the lowest QNM of the black hole. After about 
t = 60M, the hole gradually settles down to its final, 
spherical configuration. 
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FIG. 7. We show the time evolution history of the Gaussian 
curvature for the highly distorted hole (Qo = 1.0). We plot 
the curvature as a gray scale using the embedded z value 
as the y— axis, and t/M as the x— axis. We note that even 
though the hole is initially very distorted with large curvature, 
it settles down to a damped oscillatory pattern at later times 
with a frequency of about 17M. 

In Figs. |^a and |^b we show a similar diagram for 
the two black hole collision (Misner /x = 2.2) evolved in 
"Class I" coordinates. At late times, k has an appear- 
ance similar to that shown by the highly distorted case 
discussed above. At early times, there are two separate 
black holes whose horizons are about to merge. The sur- 
faces are most highly distorted along the caustic line, 
and the Gaussian curvature is largest there. In Fig. ||a 
we show the entire history of k for this system. We see 
that in the early times, before coalescence, the Gaussian 
curvature is very high near the coalescence point, (k is 
in fact singular on the EH at the caustic; we show the 
curvature very close to the caustic). In Fig. ||b we show 
the early time behavior, so the details of the curvature 
can be seen around the coalescence point. 
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FIG. 8. We show the time evolution history of the Gaussian 
curvature for the 2BH collision (Misner /j, — 2.2). In Fig. 
(a) we show the entire history of the horizon curvature, and 
note the repeated oscillation pattern, as in Fig.[j], here with a 
more complicated pattern, but still with a frequency of about 
17M. In Fig. (b) we show the early time behavior of the 
system, seeing the strong curvature near the systems cusp as 
the holes come together. We note that the line of caustics is 
the region inside the two holes before coalescence, as indicated 
by the arrows. 



4- Embedding Diagrams and Embedding Histories 

The use of embedding diagrams to study the intrinsic 
geometry of spacetimes is not new in relativity. It is a 
particularly useful way to study a curved 2D surface on a 
constant time slice. The embedding technique creates a 
fictitious 2D surface in a flat 3D Euclidean space with the 
same geometric properties as the original 2D surface in 
the curved 3D space. This technique has been described 
fully in Ref. where it was used to study AH surfaces. 
We follow the same embedding approach here. 

We can perform a non-trivial test of our embedding 
treatment by embedding the analytic Kerr horizon, and 
comparing this horizon with the embedding diagrams 
predicted in Ref. [^7| . For high-rotation Kerr black holes 
(a/M > \/3/2), it is not possible to embed the horizon 
around the pole. In Fig. |9| we show how our EH finder 
finds and embeds the correct analytic horizon up to the 
9 value where the embedding no longer exists. 
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FIG. 9. We compare the embedding of the Kerr horizon 
from our horizon finder with the known embedding for the 
value a/M = 0.877, at which the entire horizon cannot be 
embedded into flat space. We notice the agreement between 
our finder, the diamonds, and the known solution, the line. 



In Fig. |10ja-|10|c we show a time sequence of the em- 
bedding diagram for the black hole studied in Fig. ^. In 
Fig. |lO| d we show the EH and Schwarzschild embeddings, 
where both embeddings are normalized by the horizon 
mass. This normalization removes the spurious area 
growth caused by errors in the numerical spacetime, as 
seen in Fig. ||, from our horizon embeddings. In Fig. [To| a 
we see that the initial embedding is very prolate, in con- 
cordance with the large value of C r shown in Fig. |^ at 
t = 0. We also note that the final state is indeed a 
Schwarzschild-like horizon, namely a spherical black hole 
characterized solely by its mass. 
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FIG. 10. We show a sequence of embedding diagrams for an 
EH distorted by a large amplitude Brill wave. Even though 
the horizon geometry is very non-spherical at t = 0M, as 
demonstrated by the cigar shaped event horizon, the system 
quickly becomes fairly spherical, as shown in the time snap- 
shots in figures (a), (b), and (c), and also in Fig.[|, which 
corresponds to this system. In Fig. (d) we see that the late 
time horizon is essentially a Schwarzschild horizon, as it has 
settled down to a sphere with radius 2M. 



In Figs. |1 l|a-|l l|d we show a time sequence of embed- 
ding diagrams of the EH for a Bowen and York, rotat- 
ing black hole 32 2j|, with angular momentum J = 15, 
evolved by a code described in Ref. ||. We see from 
Fig. [ll]a that the initial EH is quite spherical (We make 
the front 45° of the horizon transparent to facilitate view- 
ing, hence the "pacman" appearance of the horizon) . The 
Bowen and York construction differs from that of a sta- 
tionary Kerr hole, so this data set can be regarded as 
containing a gravitational wave that makes the initial 
black hole horizon more spherical than the oblate pure 
Kerr hole. At t = 12AM into the evolution, as shown 
in Fig. p"l|b, the embedding has a shape reminiscent of a 
napkin holder; the top and bottom sections of the hori- 
zon cannot be embedded as they have negative curvature 
at the axisymmetric pole, and therefore cannot be repre- 
sented in a Euclidean space. At this instant in time, the 
extent of the unembeddable region is near a maximum. 
Fig. [TT| c shows the geometry at a later time, as the hori- 
zon settles down towards its Kerr form. Fig. [TT| d shows 
a quadrant of the EH embedding at time t — 45M . At 
this time, the hole has settled down to the Kerr form, in 
accordance with the no hair theorem ^T|. The shape of 
the EH is to high accuracy the same as that of an ana- 
lytic Kerr hole of a/m — 0.877, the embedding of which 
is plotted as a dotted line for comparison with the nu- 
merical result. We note that the value of a/m = 0.877 
is the value of the rotation specified in the initial data 



solve (a J = 15 Bowen and York hole), and that this 
result is still observed late in the evolution. This is phys- 
ically required, as an axisymmetric system cannot radi- 
ate angular momentum. The fact that our horizon finder 
confirms this late time behavior is a strong verification 
of the accuracy of our methods. Notice that there is still 
a region of the horizon that cannot be embedded, as the 
horizon for such rapidly rotating black hole is "too flat" 
for Euclidean space, and the regime in which the EH can- 
not be embedded matches the region for a Kerr EH, as 
was also shown in Fig. 0. 
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FIG. 11. We show a sequence of embedding diagrams of the 
EH for a Bowen and York rotating black hole, with angular 
momentum J = 15 (a/M = 0.877). We show snapshots of 
the horizon at (a) time t = (b) time t = 12.4, and (c) time 
t = 34.5. To allow clearer visualization of the region where 
the embedding fails, we make the front 45° of the horizon 
transparent, hence the "pacman" appearance of the holes. We 
note that, at late times, the system approaches the analytic 
Kerr embedding. 



In Figs. |12|a— |12|d we show four snapshots of the em- 
bedding of the EH for the two black hole head-on col- 
lision case generated with the Cadez coordinate system 
for /i = 2.2. Fig. 12 a shows the embedding of the EH 
on the initial, time symmetric slice (t = 0). We see 
the two individual black holes, with cusps on each hori- 
zon on the z— axis. In Fig. 12 3 we show the embedding 
at time t = 5.4, shortly after the merging of the two 
holes. In Fig. |l2| c we see the late time spherical behavior 
of the horizon, despite the system's tumultuous begin- 
nings. In Fig. |l2|d, we compare the embedding of the EH 
at t = 80M, shown as a solid line, to the horizon of a 
Schwarzschild hole with the appropriate mass, shown as 
a dashed line (once again, we normalize our final embed- 
ding by the final area). Again we see the no hair theorem 
at work, in that the initial condition with no charge and 
no angular momentum settled down to a black hole com- 
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FIG. 12. We show a sequence of embedding diagrams of 
the EH for the two black hole head-on collision case at (a) 
time t = (b) time t = 5AM, and (c) time t = 40.0A/. The 
simulation here was generated with the Cadez code. We note 
that at late times, as shown in (d), the system approaches the 
appropriate Schwarzschild hole, as expected. 

To bring out the dynamics of the horizon evolution, it 
is useful to show the "embedding history diagram" of the 
horizon instead of a series of snapshots. In Fig. |l3|, we 
show the evolution of the embedding in time for the two 
black hole case just discussed. In this diagram the 4> di- 
rection has been suppressed, i.e., we stack up <p = const 
cross sections of the 2D embeddings from various times 
to create a continuous, 2D embedding history diagram. 
We note that this figure is not a spacetime diagram, in 
that the (p, z) space away from the horizon surface has no 
physical or mathematical connection to the curved 3 + 1 
spacetime. However, these embedding history diagrams 
are a convenient and effective method for showing the 
evolution of the embedding of the event horizon surface 
in coordinate time (t) in the fictitious Euclidian (p, z) 
space. This figure shows the geometry of the individ- 
ual holes as they approach each other, with a cusp on 
each horizon. The distance between the holes before the 
merger, which is not prescribed in the embedding process 
for data generated in the Cadez coordinates, is chosen to 
keep the embedding history diagram smooth. After the 
merger, one can (barely) see the oscillation of the final 
horizon, which occurs at the normal mode frequency of 
the final black hole. In this diagram we also show the 
evolution of various horizon generators (photons moving 
normal and tangent to the horizon) as lines on the sur- 
face. The determination and use of these generators will 
be discussed in detail in the next section. 




FIG. 13. We show the embedding history diagram of the 
two black hole collision generated with the Cadez code. This 
diagram, the famous "pair of pants" diagram, shows a time 
history of the embedding of the horizon by stacking consecu- 
tive embeddings on top of each other in time. The lines on the 
surface show the paths of the horizon generators, and shows 
them leaving the surface at the crossover caustic, as will be 
discussed below. 

Another interesting embedding history diagram is 
shown in Fig. [l4|. Here we show the embedding of the 
equator of the horizon of a Bowen and York black hole 
from Fig. [ll]. We see the equator bulge out and then back 
in, as the hole becomes more and less prolate (the total 
area increases in time). As discussed in Sec. IV we em- 



bed the equator since it allows us to show the generator 
motion in the <b direction. 



10 




FIG. 14. We show the embedding history diagram for the 
case of a distorted Kerr black hole. Here we suppress the 6 di- 
rection, and embed only the equator at all <j> values. Although 
the spacetime is axisymmetric, so there is no (j> variance, this 
embedding demonstrates the rotation of the generators in the 
<j> direction as the system evolves in time. This diagram is a 
numerical construction of the "barber pole twist" diagram. 



C. Numerical Convergence of the Horizon Measures 

The study of numerical convergence is important for 
any numerical treatments based on finite differencing, 
and we discuss it for each of our results here. We give 
a brief overview of numerical convergence here, but for 
a more detailed discussion, see Refs. [ p3|]34| ]. To avoid 
confusion with Paper I, we emphasize that the numerical 
convergence we discuss here is the usual convergence rate 
of our numerical results depending on grid resolution. It 
is a completely different phenomenon than the physical 
convergence discussed in Sec. || and in Paper I, which 
is a physical attraction of null surfaces to the horizon 
independent of numerical treatment. 

Given three solutions to a discretized equation, L, M, 
and H at resolutions Ax, Ax/q, and Ax/q 2 , the conver- 
gence exponent is defined as 



where — is simple subtraction for numbers, and a combi- 
nation of interpolation onto a common grid and reduction 
via a norm operator for fields. The measure a indicates 
that the error in a numerical solution is of order Ax° . 

In Fig. [15] we show the numerical convergence expo- 
nents of the horizon area A and ratio of circumferences 
C r in a slightly distorted single black hole evolving in 



time (Qq = 0.1, r] = 0, Co = 1.0, n = 2). We choose 
this case for our convergence studies since the spacetime 
is quite accurate so effects of numerical error in the back- 
ground spacetime are minimized, and we can directly test 
our horizon treatments (We see similar convergence re- 
sults for all the spacetimes discussed here; since we have 
not assumed the Einstein Equations hold in any of our 
analysis so far, constraint violations in the spacetime will 
not affect the convergence of the system, although they 
could in principle cause the horizon analysis to converge 
to a non-physical result). The convergence study is made 
by keeping the spacetime resolution fixed in all runs and 
adjusting only the number of points which represent the 
horizon. We use an interpolator of order equal to or 
higher than our evolution method on a numerical grid of 
data. 

We see that the measures A and C r converge at second 
order as expected. These quantities are simply measures 
of the interpolated metric and the surface (evolved with 
a second order MacCormack method) so any result below 
second order would signify an error. 

Additionally, we measure (but do not show) the av- 
erage convergence of the z-coordinate of the embedding 
over the entire surface (In the embedding procedure, only 
the z-coordinate is integrated; The p-coordinate is ex- 
actly given as a function of z and the metric). The em- 
bedding converges at first order. This is to be expected, 
since we use a first order integration over the derivative 
of the surface to form the embeddings. Since the em- 
bedding is only measured, not evolved, this first order 
nature is satisfactory. We note that using a higher order 
integration scheme would not substantially improve the 
accuracy of our embeddings, since we cannot remove inte- 
grals over derivatives of the surface from our embedding 
procedure. 
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FIG. 15. We show the time evolution of the convergence 
exponent of A and C r in the low amplitude gravitational wave 
plus black hole spacetime. We note second order convergence 
throughout the entire run, and that the convergence rate of 
C r has (small) spikes associated with an oscillatory function, 
while that of A does not. 
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Notice that with this choice, the null geodesic is not 
affinely parameterized, but instead, adapted to the global 
time coordinate t used in the numerical calculation of the 
spacetime itself. 

One important advantage of determining the null gen- 
erator using Eqs. ([t2]) and ( |l3| ) is that in this formulation, 
the trajectories obtained are guaranteed to lie on the EH. 
This is in contrast to numerically integrating the second 
order geodesic equation directly. As shown in Paper I, 
integration of the geodesic equation directly can lead to 
spurious tangential drifting effects which can significantly 
affect the position of the horizon generators. This differ- 
ence can lead to errors of interpretation, as described in 
Paper I. The importance of obtaining accurate trajec- 
tories of the horizon generators is clear. Generators of 
the horizon contain all the information of the dynamics 
of the EH. The entire membrane formulation described 
Sec. [y] is based on these trajectories. Thus, inaccurate 
location of the generators due to tangential drifting can 
make analysis of the horizon dynamics via the generators 
impossible. 



B. Analytic Kerr as a test case 



IV. HORIZON GENERATORS 

We have already seen examples of generators of the 
horizon in several of the above figures. In this section we 
show that the horizon generators can be located in a nu- 
merical spacetime using information already constructed 
in our surface-based horizon finder. We will use these 
generators to study the motion and dynamics of black 
hole event horizons. 



A. Formulation 

The EH is generated by null geodesies. With the EH 
given by f(t, x r ) = 0, the null geodesies that generate the 
surface satisfy 



dx a 
~df 



(12) 



where A(x^) is a scalar function of the four coordinates. 
Notice that in terms of /, the generators satisfy a first 
order equation, rather than the more complicated second 
order geodesic equation. We choose the normalization 
A{x> i ) to be 



A{x») 



(13) 



so that the null vector tangent to the null geodesies is 
given by 



We briefly study the motion of generators in the an- 
alytic Kerr case. In this generator will rotate in 
the (^-direction on the horizon with a rate 



d<t> 
~dt 



2M 2 + 2M{M 2 



2)1/2 



(15) 



In Fig.[l6| we show d<j>/dt for various Kerr spacetimes. 
We measure the <j> location of the horizon generators, 
numerically differentiate in time, and show the result as 
solid black diamonds. We compare these results with the 
analytic result, Eq. [l5|, shown as a solid line. We note 
the results agree with the analytically expected value. 
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FIG. 16. We show horizon generator angular velocity, 
d<f>/dt for a M — 1 black hole with values of a/M between 
and 1 as solid diamonds. We compare these results with 
the analytic value, shown as a solid line. We note the excel- 
lent agreement 



C. Horizon Generators in Dynamical Spacetimes 

We now apply these techniques to the study of the tra- 
jectories of the horizon generators for three numerically 
constructed dynamical spacetimes. 

We first consider the low amplitude Brill wave plus 
black hole spacetime considered above (Qo = 0.1). In 
this spacetime we expect a non-spherical evolution in 
the generators. Rather than just moving radially, as the 
generators would in a dynamically sliced spherical space- 
time, we also expect some non spherical deflection to be 
noticeable in the generators. We show this deflection 
by plotting the difference between the generator angu- 
lar location, 9 gen and the late time generator location, 
versus 9q itself, evolving in time. Equatorial plane 
symmetry requires there is no deflection at the equator, 
and axisymmetry require that there is no deflection at 
the pole, thus all the generator deflection must occur 
between the equator and the pole. In the intermediate 
region, the generators oscillate with a quasi-normal mode 
frequency with an amplitude dying down at late times. In 
Fig.[l7], we show the deflection quantity 9 gen — 6q evolving 
in time, and note that the deflection is small, but displays 
this expected behavior. 




FIG. 17. We show the angular deflection of generators for 
a horizon with a low amplitude brill wave initially incident on 
a black hole. We show the deflection by plotting the angular 
location of the generators, (9 gcn minus their late time position, 
8o, versus their late time position, 9o, evolving in time for all 
generators. From the figure it is clear that the angular de- 
flection occurs away from the equator and pole, as is obvious 
from simple symmetry arguments. 

In Fig. H, the "Barber Pole Twist" diagram, we have 
seen the <fi motion of the generators in Kerr-like space- 
times. In Figjl^, we plot the quantity dip/dt of the pho- 
tons versus time. We see that they settle down to a 
constant value at late times, with 



dt 



(t = BOM) = 0.293. 



(16) 



This is to be compared to the analytic value of 0.296 given 
by Eq. @ with a/M = 0.877, denoted by the dashed 
line in the figure. We see that the measured value at t — 
50M differs from the analytic value by about 1%, which 
demonstrates that the hole is settling down to a Kerr 
black hole at late times, and that the determination of the 
horizon generators is quite accurate, although late time 
errors in the numerical spacetime lead to the observed 
small difference from the expected value in the generator 
angular velocity. 
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FIG. 18. We plot the generator angular velocity d<j>/dt for 
a horizon generator on the equator vs. t for the distorted 
Kerr black hole. We note that, although this quantity is not 
constant in time, it approaches the analytic value of 0.296 at 
late times, as the horizon settles down to its Kerr form. 

Turning to the "Pair of Pants" diagram, Fig. [l^, which 
shows the embedding of two colliding BHs, we see that 
the most interesting feature of the generators is that some 
of them leave the horizon (going backwards in time) at 
the inner seam of the pants. There is a line of caus- 
tic points on the z-axis extending backward from the 
"crotch" point where the two horizons merge. It is at 
these points along the caustic line in the history diagram 
that photons originally travelling in the causal past of 
null infinity (J~(T + )) join the horizon as generators. As 
discussed above, only the surface of the horizon has been 
embedded; the photons that have left the embedding dia- 
gram have also left the embedding space, and their paths 
are only shown to denote their joining the horizon. 

In Fig. 19 we show the coordinate location of the gen- 
erators and horizon surface found using the Cadez code. 
The EH location at various times is shown by heavy solid 
lines. The t — surface is the horizon of two distinct BHs 
at the initial time, which evolves to a single, merged hori- 
zon, shown at t = 3.1M. We see that generators which 
start outside the EH (denoted by inward pointing trian- 
gles in the figure) move inwards, cross on the z-axis, and 
join the horizon. This crossing of generators of the EH 
in the two black hole collision is crucial to recent under- 
standing of the structure of the horizon in the Misner 
spacetime. Further analysis of the nature of such lines 
of caustics is possible and underway pq ]. Coupled with 
new techniques for evolving multiple black hole space- 
times, our techniques should allow an increase in our un- 
derstanding of how the generators behave in dynamical 
multiple black hole spacetimes. 
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FIG. 19. In this diagram we analyze the trajectories of the 
photons before and after they join the horizon for the case 
of two colliding black holes (Misner /i = 2.2). Slices of the 
horizon are shown at t = 0, 1.9, 2.7, and 3.1M. We note that 
generators not originally on the horizon (shown by inward 
pointing arrows on this figure) cross over each other at a line 
of caustics on the z— axis and join the horizon as the holes 
collide. For example, the two photons labeled A and A' join 
the horizon at t — 2.7 M, crossing over at the point shown as 
an open circle. At time t = 3.1M they are on the horizon at 
the points shown as solid circles. The second black hole is not 
shown, as the system has equatorial plane symmetry. 



D. Numerical Convergence of the Generators 

We can measure convergence of the generator loca- 
tions, just as we measured convergence of horizon mea- 
sures. Since the generator location is an ODE integration 
with coefficients determined by the surface location and 
the derivatives of the surface, the appropriate test is to 
keep the number of generators fixed, while changing the 
spacing of the surface. We can then measure the differ- 
ences in generator locations as a function of spacing of 
the surface and form a convergence measure for each gen- 
erator, which can then be averaged over all generators. 

We show the result of performing this operation on the 
radial and angular positions of the generators in Fig. |2(], 
using the low amplitude Brill wave plus black hole space- 
time considered above. We note that the radial position 
of the generators (solid line), which is non-oscillatory, 
converges at second order. However, the angular position 
(dashed line) has spikes typical of an oscillatory func- 
tion, but converges below second order. This lower order 
convergence is due to the principal term in the angu- 
lar generator position evolution being the (interpolated) 
derivative of the horizon surface. That is, since we intcr- 
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polate second order spatial derivatives of the surface for 
the generator sources, the evolution of the generator an- 
gular positions has error terms larger than the Aq 2 terms. 
This convergence order could possibly be increased by us- 
ing fourth order spatial derivatives and very high order 
interpolators. 
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FIG. 20. We show the convergence of r and 6 for the hori- 
zon generators in the low amplitude Brill wave spacetime. 
The convergence exponents a r and ag are plotted versus time. 
We note that the radial location of the generators converges 
at second order, but the angular location converges at some- 
what less than second order. This is unsurprising since the 
evolution equation for generators' angular location is domi- 
nated by an interpolation of the numerical derivative of the 
horizon surface. 



V. THE MEMBRANE PARADIGM 

We now turn to a detailed analysis of the information 
carried by the congruence of horizon generators, and the 
extent to which this can be used in numerical relativity 
as a tool to investigate black hole dynamics. The the- 
oretical basis for this study is based on the membrane 
paradigm (MP) g]. The MP views the black hole as a 
2-surface in a 3-space with the properties of a viscous 
fluid. In many ways, the EH in a dynamical spacetime 
is like a soap bubble perturbed by external influences. 
The MP is particularly valuable in providing an intuitive 
understanding of how a BH reacts to its surroundings. 

There has been much study of gravitational interac- 
tions using the MP in quasi-stationary situations |56],[37j . 
With the advent of numerical identification of the EH 
and generators described above, we can now start to con- 
sider applying the MP to fully non-linear and dynamical 
spacetimes. With this goal in mind, we demonstrate how 
to construct the MP quantities on a numerically located 



EH, and examine the accuracy of these constructions in 
several testbed spacetimes. 



A. Formulation 

We begin by discussing the MP formalism with the goal 
of being able to construct MP quantities on our numeri- 
cally located horizons. The membrane paradigm requires 
the choice of a time slicing, splitting up spacetime into 
an "absolute space" and a "universal time" JjJ . To apply 
the MP to numerical relativity, we choose the universal 
time to be the same as the time coordinate t used in 
the numerical evolution. This implies that the time co- 
ordinate used in the numerical evolution has to be well 
behaved on the EH. This is the case for all of the black 
hole spacetimes we have numerically constructed. 

We define the four vector I to be the tangent to the 
horizon generators, and we normalize it as in Eq. ( |T4| ) 
above, with t being considered as the "universal time" . 
This vector is in the full 4-dimensional space, which we 
index with Greek letters, = (0,1,2,3). On the 

2D spacelike section of the EH at constant t, we choose 
spacelike 2D coordinates x a which we index with lower 
case Roman letters, a, b, ... = (2, 3), which are comoving 
with the horizon generators, i.e., 



d_ 

dt 



d_ 

dt 



(17) 



where i is the comoving generator time coordinate (which 
is identical to the time in the simulation by Eq. (|14|)). In 
a coordinate basis, we have the spatial basis vectors 



d 

dx a 



(18) 



which are orthogonal to £ by construction. We define the 
fourth basis vector ft by 



ft ■ ft = 
ft-t = -1 
ft ■ e a = 0. 

The induced metric on the 2D horizon section is 

lab — e a e b9nv 



(19) 



(20) 



In the membrane paradigm the description of the dy- 
namics of the horizon is given in terms of the horizon 
surface gravity gn, the shear a^ b , the expansion , and 
the Hajicek field Q„ . They are defined as 



ff =iiLlndet7" 
2dt ' 



rr»- 1 

°ab-^ 



dt lab 



nf = -ft ■ v a i 



(21) 

(22) 

(23) 
(24) 
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These quantities are dependent on the choice of time 
coordinate t, as they explicitly involve t in their defini- 
tion. That is, they are gauge dependent measures of the 
horizon dynamics. In the formulation of the membrane 
paradigm given in Ref. [|l| , a particular time slicing is cho- 
sen for a stationary black hole, e.g., a Kerr black hole. 
In this slicing, without perturbation, gn and Qh take 
on special values while Q H and a^ b vanish. For small 
perturbations about a Kerr horizon, gn and are first 
order slicing dependent. In the formulation given in Ref. 
0, time slicings of the perturbed black hole are cho- 
sen so that the surface gravity gn remains unchanged in 
time. In our application of the membrane paradigm to 
numerical relativity, as we are mostly interested in highly 
dynamical and fully nonlinear interactions, we do not put 
such restrictions on the time slicing. Rather, we let the 
time slicing be determined by the natural choice of the 
numerical evolution (maximal slicing for most cases pre- 
sented in this paper). We expect that the new features 
introduced by different slicings will become familiar when 
the formulation is used in more black hole studies, and 
hopefully allow further insight into the slicing conditions 
and numerical evolutions. 

The horizon quantities ( pl| ) - ( |2~4| ) satisfy the following: 
The "tidal force equation" 



the "focusing equation" 



-£. 



H 



(25) 



1 



Dt&H = 9hQ H ~ - Vab^H - 8tt V-T, (26) 

and the "Hajicek equation" 

An? + (of c + ~<^e H )nf + erf = 

{gh + \Q H ),a-^\\ b + ^T all t. (27) 



A =_L I ■ V is the projection of the covariant derivative 
along £ into the horizon section. "||" denotes covariant 
differentiation on the horizon section. C M „ pcr is the Weyl 
tensor and T^ v is the energy-momentum tensor. 

The comparison of Eqs. (|25|)-(|2^) with the evolution 
equations for a 2D viscous fluid gives meaning to the hori- 
zon quantities Eq. (§l|) - @). One finds that Eq. (§|) 
describes the response of a fluid to a gravitational tidal 
field, Eq. (|2^) describes the energy conservation of the 
viscous flow, and Eq. ( p7| ) is the corresponding Navier- 
Stokes equation of the fluid flow. The surface density of 
the mass-energy of the fluid is identified as — ©#/87r, the 
surface pressure is gu/^, and the momentum density 
corresponds to — il^ / 87r. The dynamics of the EH of a 
black hole can be understood in analogy to the motion 
of a fluid on a soap bubble. In the following section, we 
show how these "fluid" quantities can be constructed for 
an EH located in a numerical simulation. 



B. Constructing Membrane Quantities 

Once f(t,x l ) = is given, we obtain £ as given in 
Eq. (|l4]) in a straightforward manner. Next, we define 
the comoving coordinates x a , (a — 1,2) on the horizon 
section by (9, cf>) . Then we have 



e a = (d§, d$) =(j>,q). 



(28) 



The coordinate components of the two basis vectors can 
be obtained by 



p = d g =p r d r + P e d 9 +p%, 

where p r is defined to be 

r dr 
P = d§ 

difference in r for neighboring generators 
difference in 9 for neighboring generators ' 



(29) 



(30) 



and likewise for p , p^. We use this definition in a dis- 
crete fashion, differencing over generator locations, and 
therefore our basis vectors will always have a discretiza- 
tion error based on the initial spacing of generators in 
9 space. As the coordinates 9 and 4> are chosen to be 
comoving, we have = = (f. For the axisymmetric 
cases considered here, we pick <fi = <p and thus q = d/dcf>, 
the azimuthal killing vector. 

The horizon two -metric is then written as 



lab 



i fin 



7 « 

' ftfh 



7, 



H 



H 

e<t> ><p<p 



The individual components are defined by, e.g., 



7l = gijP'P 3 



(31) 



(32) 



Solving for fi is particularly troublesome. We use the 
following, geometrically motivated, method. When we 
solved for d t f in Eq. ||), we solved the quadratic equation 
choosing the positive root for the outgoing null surface. 
We could also have chosen the negative root, and found 
an evolution equation for the ingoing null surface. Let us 
call the ingoing evolution equation <9 t - /, and Eq. (^) <9 t + / 
temporarily. We will use the notation d^f — (df f,dif). 

Thus, we can form two null vectors L and N as 



L lt = (8tf,d i f) 

N li = (drf,d i f). 



(33) 



From Eq. (g) it is clear that both L and N are null, and 
that L is simply I with a different normalization. 

However it is also clear that N ■ e a = 0. To see this, 
recall that e a has only spatial components, so 



N^ a = e^N, = eldj = e % a L t = L M < = 



(34) 
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since L is proportional to £ which is orthogonal to e a by 
construction. 

So now all that remains is to find a normalization such 
that n ■ £ = — 1. This is straightforward. Since £ = 



L/A(x u ), using Eq. (13) it is clear that 



N ■£ 



g ta d+f 



B{x^) 



and so we can define n by rescaling N by B(x^), 
n = -N/B{x^). 



(35) 



(36) 



We note we can use Eq. ( |19[ ) to measure how accurately 
n and £'s orthogonality with e a is maintained. 

Once the horizon 2-metric 7^ and full set of comoving 

vectors, (l,n,p,q) are obtained, we can form the expan- 
sion, shear, and Hajicek Field via Eq. ( pl| ) - (p4f). From 
Eq. (^i|), the surface gravity is 



(37) 



for our particular parameterization of £. 

There are several terms in the definitions of the mem- 
brane quantities which require careful numerical and 
analytical treatment in order to be evaluated in our 
framework. In particular, in order to evaluate the hori- 
zon quantities accurately, we must be able to evalu- 
ate d"/ a b/di, preferably without taking numerical time 
derivatives. From Eq. (|32"|), the horizon 2-metric 7 has 
two types of terms, those due to the comoving basis vec- 
tors p and q, and those due to the spacetime 4-metric, 
cjij. Thus using the chain rule to evaluate djab/di will 
yield terms like dgij/dt and dp z /dt. 

The derivatives along the generators of the spacetime 
metric can be evaluated using the metric evolution equa- 
tions. We note that this is the first point we have used the 
evolution equations for </y, and therefore the accuracy 
with which our spacetime obeys these evolution equa- 
tions enters can enter into our quantities. In other words, 
if the relationship between dtg and 2aK is only obeyed 
to a given order, we cannot expect our quantities which 
use this relationship to be obeyed at a higher order. By 
virtue of £* = 1, 



d 9ii = , Mg 



9ij 



-2aK {j + Dfpj + Djfa + £ k g ij<k (38) 



where Kij is the extrinsic curvature of the 3 surface. 
and g-ij are both readily available in the numerically con- 
structed spacetime. 

We can find the terms dp 1 jdt by commuting partial 
derivatives. Namely, 



dp 1 
~W 



d dx l _ d£ l 
I¥ ~ ~d9 



(39) 



The time derivative of p is the spatial derivative of L We 
can evaluate the spatial derivative of £ with a single time 



slice finite difference of our surface and surface quantities, 
and thus find the required time derivatives. To summa- 
rize, expanding the time derivative of the horizon metric 
using the chain rule, and using the above two techniques, 
we can find the d^ a b/dt terms in a single time slice. 

Thus, we have a method for finding the four horizon 
quantities which describe the kinematics of the horizon 
surface. This method is contained entirely in a single 
3-slice. We should note that it is also possible to cre- 
ate the membrane quantities in a direct fashion using 
numerical derivatives in time to evaluate the expansion 
and shear. We call this the "time difference" evaluation 
of the expansion, as opposed to the "single slice" eval- 
uation. We find that the single slice method invariably 
gives smoother and more accurate data for the membrane 
quantities than the time difference method. 

An additional difficulty comes in evaluating the hori- 
zon equations, Eq. (J2q)-(|27|). Two terms pose a difficulty 
there, D^Z a and &a\\bi where Z is any tensor on the 
horizon. Luckily, we only need to evaluate these terms 
as a check; we do not use the horizon equations in our 
evolution. Thus we can use first order accurate methods 
to evaluate these if need be. 

We first turn our attention to D^Z a . First we introduce 
a Christoffel symbol for the (t, 9, (j>) coordinates (e.g., the 
null horizon 3-surface in co-moving coordinates, which 
we will here index with (q, r, ...)). We denote this as 
( 3 )F« ,. We find 



Di z a = d -^-^ ai z q . 



(40) 



The horizon "3-metric", 7^ is simply given by j a b if 
q,r 7^ t and elsewhere. Thus ^ 3 ^r 9 rs can be simply 
evaluated as 



Thus we can easily evaluate Eq. (Eot) as 



DfZa 



dZ a 
St 



1 



l bC lba,tZ a . 



(41) 



(42) 



The only term which we cannot calculate in a single slice 
is dZ a /di, but we can simply calculate that by storing the 
quantity Z a at three time steps and then use a centered 
time derivative to evaluate the term at the middle step 
after all three steps are taken. 



The term <J a \\b is evaluated directly, e.g. 



„ b 
Va \\b 



a a h 



(2)p6 



(43) 



and the 4 independent non-zero terms of ^T a b c are eval- 
uated directly from spatial derivatives of the horizon 2- 
metric. 



17 



C. Test and Applications of MP Quantities 

In this section we apply the membrane quantities to a 
set of testbed analytical and numerical black hole space- 
times that have been computed using codes described in 
Refs. Hl^]. Our aim here is to probe whether these tools 
can be used in a practical way to explore the dynamics of 
black hole horizons in numerically generated spacetimes. 
We will consider the physics of these quantities, for a set 
of interesting spacetimes, in a future paper. 



1. Flat space 



Flat space in Minkowski coordinates, 
ds 2 = -dt 2 + dr 2 + r 2 dn 2 



(44) 



allows us to test our expressions for against easily un- 
derstandable analytic solutions. Although flat space has 
no EH, it does have null surfaces, and our construction 
carries over to them. 

Most notably, we know that for spherical null surfaces 
in flat space, the expansion of a sphere of radius r is 



A dt r 



(45) 



since in flat space, A = Anr 2 and dr/dt = c = 1. Using 
this relationship, we can trivially check our expressions 
for 8. Additionally, we can form dA/dt from integrals 
of the expansion, which carries over into the dynamical 
black hole case, where we can compare this integral of 
<d H with a numerically calculated dA/dt. Evaluating 
the expansion in flat space gives the expected answer. 



2. Analytic Schwarzschild 

We next turn to the analytic Schwarzschild spacetime 
described in standard coordinates, 



ds 2 



2M 



dt 1 



dr 2 



+ r 2 dn 2 . (46) 



In this spacetime, the expected results are that the gen- 
erators and surface will be attracted backwards in time 
towards the true horizon (at r = 2M), that Qh, &ab 
and approach zero exponentially as the surface ap- 
proaches the true horizon, and that gn approaches the 
analytic value of 1/2M. Moreover, we can check that re- 
lationship between the integral of the expansion and the 
area change holds in this spacetime by numerically dif- 
ferentiating the horizon area, which allows another test 
of our expressions. These relationships are obeyed. 

In analytic Schwarzschild we can trivially evaluate the 
above expressions for £, and gu on an arbitrary null 
sphere of radius r to find 



2M 

= 1,1 ,0,0 



r r \ r 



2M lH r 



2M 



9H 



1 



2M 



(47) 

(48) 

(49) 
(50) 



Note that l r and vanish on the horizon (r = 2M) as 
expected, and gu takes the value 1/2M . We check these 
relationships for surfaces away from the horizon and we 
see that our surfaces give the analytic results for all null 
spheres in the spacetime. 

Additionally, each of the horizon equations, Eq. (p5|)- 
( p7| ) should be obeyed in this spacetime. We evaluate 
only the focusing equation violation, however, since the 
tidal force equation contains the electric part of the Weyl 
tensor, £K, which causes this equation not to be a check 
on the membrane quantities alone, and the Hajicek equa- 
tion is trivially satisfied with a spherically symmetric gu 
and Of = 0. 

The vanishing of the focusing equation violation allows 
us a strong check on our method. Since the focusing 
equation requires the covariant derivative of the expan- 
sion, D t 0, we expect the focusing equation to be obeyed 
as accurately as D t Q is evaluated. Recall, we evaluate 
D t <d by taking a centered finite difference in time, so 
we expect the focusing equation violation in our space- 
time to converge towards zero at 0(At 2 ). We test this 
by finding a surface in the analytic Schwarzschild back- 
ground first using a Courant factor A = 0.2 and then 
A = 0.4, doubling the time step. We then measure the 
focusing equation violation in these two runs. If the re- 
sult is converging towards zero, the focusing equation 
violation should be four times larger in the A = 0.4 case. 
We demonstrate this convergence in Fig. ^l] by plotting 
the focusing equation violation with A = 0.2 as a line, 
and by plotting one quarter the focusing equation viola- 
tion with A = 0.4 as diamonds. The demonstration that 
these two sets of data are the same indicates that we are 
converging towards a surface which satisfies the focusing 
equations. We note that, as the surface becomes very 
close to the actual horizon, the focusing equation is zero 
at levels close to machine precision in both simulations, 
so convergence can no longer be observed numerically. 
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t/M 

FIG. 21. We show the violation of the focusing equation, F, 
for a given sized time step as a solid line, and one quarter the 
violation for double the time step as diamonds, for a surface 
integrated in the analytic Schwarzschild spacetime. The fact 
that these data are coincident indicates that we are converging 
towards a null surface which satisfies the focusing equation. 
The surface in question starts at r = 2AM at t = 30. The 
exponential shrinking of the violation is directly due to the ex- 
ponential approach of the expansion towards zero. Note also 
that as the focusing equation approaches machine precision 
levels (here 10~ 14 ) convergence fails, since both quantities 
are effectively zero. 



3. Maximally Sliced Schwarzschild 



With the advent of new hyperbolic systems for the 



Einstein equations 
ary Conditions |24 



38| |40| and Apparent Horizon Bound- 
HF, long time highly accurate one di- 
mensional evolutions of a maximally sliced Schwarzschild 
black hole are quite readily available, and so we can use 
these very accurate spacetimes to test our horizon finding 
method. In this section, we consider a maximally sliced 
black hole evolved with the eigen-method code described 
in Ref. |39f| , which allows long time evolution with excep- 
tionally small error. 

We first can test the evaluation of the horizon 2-metric, 
jab- In the case of no angular generator motion, where 
the generators are chosen to be identically on the points 
on which the horizon surface f(t,x l ) is evolved, the hori- 
zon 2-metric 7 a fc and the induced surface 2-metric used 
to evaluate area and circumferences should be identi- 
cal. That is, we should get the same answer evaluat- 
ing Eq. (R) whether we use j a b as defined by Eq. (||) 
or Eq. (pfj|). Moreover, the vectors p and q should have 
components (0, 0, 1, 0) and (0, 0, 0, 1) respectively. We see 
both of these features to machine precision in the maxi- 



mally sliced Schwarzschild spacetimes. 

Spherical symmetry also leads to a vanishing shear; our 
expression for the shear vanishes to machine precision. 
However the expressions for the expansion is non-trivial, 
and since we have a very small (but non-zero) area growth 
due to numerical error, we can very accurately measure 
how well the expansion measures area change. 

We choose two trial surfaces for our test, one slightly 
outside the horizon and one slightly inside, and integrate 
them backwards in time. As expected from Paper I, these 
surfaces converge towards each other rapidly, and lock 
onto the same surface, but have some non-trivial area 
change, due to the "locking on" process before the sur- 
faces join the horizon, and due to numerical error after- 
wards. In Fig. 22 we plot dA/dt calculated by differen- 
tiating the area reported by the code, and also by inte- 
grals over the surface of the expansion. We see that these 
quantities agree, strongly indicating that our evaluation 
of the expansion is correct. 
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FIG. 22. We show dA/dt evaluated by taking both the 
numerical derivative of the area calculated by the code, and 
surface integrals of the expansion found from the comoving 
horizon two metric. We use a very accurate maximally sliced 
Schwarzschild spacetime which has a very small, but non-zero, 
numerical error in the spacetime. We integrate two surfaces, 
one originally inside and the other originally outside the event 
horizon. We note the excellent agreement between the two 
measures of dA/dt. 



4- Small Distortion Non-rotating Black Hole 

We turn to the small distortion Brill wave plus black 
hole spacetime considered above. We first test if our eval- 
uation of the horizon two metric, 7 a b, gives measures of 
the horizon geometry which are consistent with the mea- 
sures discussed in Sec. Ill B. Since the generators will ex- 
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perience angular deflection, integrals to form areas and 
circumferences will be over different coordinate locations 
when using the comoving and induced two metric. More- 
over, the measure of the geometry using the horizon two 
metric will be measured on a non-regular grid in 6, <p 
space (but a regular grid in 9, <f> space), and will therefore 
have an additional inaccuracy. Nonetheless, we see good 
agreement. In Fig.^3|, we show the difference in evalu- 
ating C p (not C r ) using the comoving and induced two 
metric in the Brill- wave plus black hole spacetime. We 
show the difference for 38 and 76 generators, respectively. 
Note as the number of generators increases (therefore re- 
ducing numerical error in the integration over the horizon 
metric due to generator deflection), the results converge 
towards the same solution, or the differences converge 
towards zero. 
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FIG. 23. We demonstrate that the generator co-moving 
metric gives accurate evaluations of the polar circumference, 
C p , in the spacetime with small amplitude Brill waves ini- 
tially on the throat. We show this by forming Cp ic from the 
induced surface metric, and C| en from the co-moving gener- 
ator horizon 2-metric. We taking the difference of the two 
measures with different numbers of generators used to form 
the horizon 2-metric. Clearly, as more generators are used 
the two methods become closer and the differences converge 
towards zero. 

We turn next to the expansion on the horizon. For the 
physical setup considered here, a gravitational wave inci- 
dent on a black hole, but with the wave centered at the 
throat, we expect the horizon to grow at t = 0, and then 
as time progresses, become static. This should show up 
as a positive expansion decreasing towards zero as time 
progresses. However, we also know that our spacetime 
has spurious area growth of the horizon due to numerical 
error in the spacetime, as found in previous studies of 
the AH. This should appear as a positive, and increas- 
ing, expansion at later times. In Fig. we show the 



expansion for this spacetime, and see exactly this be- 
havior. However, a few features of the expansion should 
be noted. Firstly, note that at late times, the expan- 
sion is not terribly smooth in time. Secondly, note that, 
near the axis (6 = 0) the expansion is somewhat oscil- 
latory. At late time and near the axis the numerically 
constructed spacetime is less accurate. We see that our 
membrane paradigm quantities as analysis tools are very 
sensitive detectors of these errors in the numerically gen- 
erated spacetime. 




FIG. 24. We show the evolution of the expansion in time 
for the horizon interacting with a small amplitude Brill wave. 
Two features are of interest here. First, we note the initial 
expansion is quite large but drops quickly, as the horizon swal- 
lows the initially incident gravitational radiation. This initial 
growth is concentrated near the equator, as the gravitational 
wave has a sin 2 form. Secondly we note that at later times 
the expansion is growing, as expected from the spurious hori- 
zon growth due to numerical error, and this growth has no 
angular dependence. We also note a small amount of noise 
on the horizon near the axis, due to spacetime inaccuracies 
there. 

This detection of error leads us to study how these 
quantities behave with changing resolution in the con- 
struction of the numerical spacetime. In Fig. ^5] we take 
the same wave parameters used above with resolutions 
of 200 x 54 and 300 x 80 to generate two spacetimes. In 
Fig. ^5], we show the area change predicted by integrat- 
ing the expansion over the 2-surface. We see that, at 
t = 0, where area change is caused by infalling gravita- 
tional radiation and the spacetime is still quite accurate, 
both systems give the same result, but at later times, 
the expansion due to spurious numerical error is consid- 
erably larger in the lower resolution spacetime, and the 
expansion appears to be converging towards zero. In the 
high resolution spacetime, the expansion is fairly inaccu- 
rate near the pole, as the system is very susceptible to 
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axis instabilities, but this noise does not show up in the 
calculation of area change, as sin 9 terms in the integral 
of the expansion kill this contribution near the pole. 
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FIG. 25. We study the behavior of the expansion with a 
moderate and high resolution numerically generated space- 
time. At early times when area growth is due to accurately 
modeled gravitational phenomena, the expansion should be 
unchanged by adjusting the spacetime. At late times, when 
area growth is due to spurious numerical error, the area 
growth should is smaller with a higher resolution simulation, 
as the spacetime is converging (at roughly second order) to- 
wards a zero-area-growth solution. 

We next turn to the shear. In this spacetime, we expect 
a non-zero shear since there is generator motion, but we 
also expect the shear <j a b to be diagonal, since the space- 
time is non-rotating and axisymmetric. In Figs. |2£] and 
p7| we plot the evolution of agg and the trace of the shear, 
a a a in time. We note that the shear is largest near the 
equator, and vanishes on the pole, as symmetry argu- 
ments require it must. (There can be no shear at the 
pole in axisymmetry, only expansion, since shear at the 
pole would imply a <j) dependence of the generator mo- 
tion). We also note that the trace of the shear vanishes 
to machine precision in Fig. ^ 
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FIG. 26. We show aee on the horizon for the low amplitude 
distortion case considered. We note that there is no shear at 
the poles, and the shear is maximal near the equator. 




FIG. 27. We show the trace of the shear, ag° + on the 
horizon for the low amplitude distortion case considered. We 
note that, even though the shear and the horizon two metric 
are order unity, this quantity effectively vanishes to machine 
precision. 



Finally, to test the surface gravity, we turn to the focus- 
ing equation, which is a complicated combination of the 
surface gravity, shear, and expansion. If this equation is 
roughly satisfied in our spacetime, then we have a strong 
verification that we are indeed measuring the membrane 
quantities appropriately. We test this by taking the aver- 
aged value of the focusing equation violation (or the LHS 
- RHS of Eq. (||^)) over the surface. In Fig. [2| we show 
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these averages evolving in time in our moderate and high 
resolution spacetimes. We note that the focusing equa- 
tion violation is small, being substantially smaller than 
the square of the shear and the expansion. However it 
is clear that the evaluation of the focusing equation vi- 
olation is also sensitive to the errors in the numerical 
spacetime and interpolations. Noise, which is generated 
from the discrete and inaccurate features of the space- 
time, is clear in Fig. |2^. However, we also observe that, 
with more spacetime resolution, the focusing equation vi- 
olation converges towards zero at approximately second 
order, as expected. 
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FIG. 28. We show the norm of the focusing equation over 
the surface in the high resolution and medium resolution 
spacetimes. We note several features. Firstly, this quantity is 
noisy, but small compared to the square of the shear and ex- 
pansion, both of which enter into the equation. Secondly, we 
note that with an increasingly accurate spacetime, the focus- 
ing equation converges towards zero at approximately second 
order. 

From the experiments in these two numerical space- 
times we conclude that our construction is appropriate 
for measuring and generating membrane paradigm type 
analysis quantities in numerical spacetimes. These quan- 
tities are sensitive detectors of error in numerical space- 
times, and they allow us to measure detailed properties 
of the event horizon and its dynamics. 



VI. CONCLUSIONS 

In this paper, we have developed a set of tools with 
which one can measure and understand the dynamics 
of event horizons in numerically generated spacetimes. 
We have shown that standard geometric measures of the 
horizon are useful tools for understanding horizon dy- 
namics. We have investigated the behavior of the gen- 



erators of the horizon in several spacetimes, including 
two black hole spacetimes, where horizons contain caus- 
tics, through which generators leave the horizon. Finally 
we presented a construction which applies the membrane 
paradigm to numerical relativity. We demonstrated that 
this construction was effective on analytic spacetimes, 
and is also applicable to numerically generated space- 
times. We also note that our techniques are applicable to 
any null surface, so could potentially be useful for study- 
ing null surface dynamics in spacetimes without black 
holes, or away from black holes. We look forward to 
more accurate dynamical black hole spacetimes, so we 
can use these quantities for detailed horizon analysis in 
numerical relativity. 
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